function b = zeroMAT(x1,x2,t,f)
%FZEROMAT Matrix version of FZERO.
% x = fzeroMAT(F,[a,b]) tries to find a zero of F(x) between
% a and b. F(a) and F(b) must have opposite signs.
% fzeroMAT returns one end point of a small subinterval of
% [a,b] where F changes sign.
% Additional arguments, fzeroMAT(F,[a,b],p1,p2,...),
% are passed on, F(x,p1,p2,..).

a = x1;
b = x2;
fa = f ( a );
fb = f ( b );
if max(sign(fa) == sign(fb))==1,
error('Function must change sign on the interval')
end
c = a;
fc = fa;
d = b - c;
e = d;
while max(abs(fb)) ~= 0
% The three current points, a, b, and c, satisfy:
% f(x) changes sign between a and b.
% abs(f(b)) <= abs(f(a)).
% c = previous b, so c might = a.
% The next point is chosen from
% Bisection point, (a+b)/2.
% Secant point determined by b and c.
% Inverse quadratic interpolation point determined
% by a, b, and c if they are distinct.
ind = (sign(fa) == sign(fb));
a   = (1-ind).*a +ind.*c;
fa  = (1-ind).*fa+ind.*fc;
d   = (1-ind).*d+ind.*(b - c);
e   = (1-ind).*e+ind.*d;

ind = ( abs(fa) < abs(fb));
c   = (1-ind).*c  + ind.*b;
b   = (1-ind).*b  + ind.*a; 
a   = (1-ind).*a  + ind.*c;
fc  = (1-ind).*fc + ind.*fb;
fb  = (1-ind).*fb + ind.*fa;
fa  = (1-ind).*fa + ind.*fc;

m = 0.5.*(a - b);
tol = 2.0*eps*max(max(max(abs(b))),1.0)+t;
if (max(max(abs(m))) <= tol) || (max(max(abs(fb))) == 0.0),
break
end
% Choose bisection or interpolation
% Bisection
indexbis = ((abs(e) < tol)+(abs(fc) <= abs(fb))>0);
% Interpolation
s = fb./fc;
indexlin = (a == c);
% Inverse quadratic interpolation
q = (1-indexlin).*fc./fa;
r = (1-indexlin).*fb./fa;
p = (1-indexlin).*(s.*(2.0.*m.*q.*(q - r) - (b - c).*(r - 1.0)))+indexlin.*2.0.*m.*s;
q = (1-indexlin).*(q - 1.0).*(r - 1.0).*(s - 1.0)+indexlin.*(1.0 - s);
ind = (p > 0);
q   = (1-ind).*q    + ind.*(-q);
p   = (1-ind).*(-p) + ind.*p;
% Is interpolated point acceptable
ind = ((2.0.*p < 3.0.*m.*q - abs(tol.*q)).*(p < abs(0.5.*e.*q)));
e = (1-indexbis).*((1-ind).*m+ind.*d)   +indexbis.*m;
d = (1-indexbis).*((1-ind).*m+ind.*p./q)+indexbis.*m;
% Next point
c = b;
fc = fb;
ind = (abs(d) > tol);
b = (1-ind).*(b - sign(b-a).*tol)+ind.*(b + d);
fb = f(b);
end